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ABSTRACT 

In this paper the capabilities of an automated CFD system which is currently available at NLR 
are demonstrated. Transonic flow around the AS28G wing/body configuration and hypersonic 
flow through a generic three-dimensional mixed-compression airbreathing inlet are simulated. An 
assessment of the level of automation of the current CFD-system is made. The problem-turnaround 
time lies within the order of a week for both applications. 

1. INTRODUCTION 

Experience with CFD technology (multiblock structured grid methods) learns that the turnaround 
time for generation of grids for Euler/Navier-Stokes calculations of complex aircraft configurations 
is large. In order to efficiently contribute to aircraft design the problem-turnaround time must be 
reduced to the order of a day or a week [1]. Unstructured grid methods offer the possibility to 
bring about this reduction. 

The level of acceptance of CFD technology in the aerodynamic design process is directly related 
to the ability to produce accurate solutions. High accuracy of aerodynamic forces is especially 
important so that the computed lift, drag and pitching moment can be relied upon to reduce the 
risks involved in aircraft design. 

In view of these aspects DLR and NLR started a cooperation entitled ”CFD for complete aircraft” 
to develop a fully automatic system for three-dimensional flow simulations. The algorithms used 
are based on the unstructured grid approach [2] which is based on a Galerkin finite-element method 
to discretise the three-dimensional Euler equations. 

The objective of this paper is to demonstrate the capabilities of the CFD system which is cur- 
rently available at NLR. The focus will be on the level of automation and accuracy of the CFD 
system. The first application concerns three-dimensional transonic inviscid flow past the AS28G- 
wing/body configuration. The second application is the study of an airbreathing inlet, which is 
part of the AEOLUS programme, a joint industry project in the Netherlands [3]. The CFD sys- 
tem is used to simulate hypersonic inviscid flow r (no real gas effects included) through a generic 
three-dimensional mixed-compression inlet with two ramps. 

2. CFD SYSTEM AT NLR 

At NLR an automated CFD-system for three-dimensional inviscid flow simulations is acquired. In 
the current CFD-system the following algorithmic steps are necessary to obtain a visual represen- 
tation of the flow field for a given aircraft configuration. 

1. Geometry definition. 

Part of this investigation has been carried out under contract awarded by the Netherlands Agency for Aerospace 
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2. Surface triangulation. 

3. Three-dimensional grid generation. 

4. Pre-processing. 

5. Flow calculation. 

6. Grid adaption. 

7. Visualisation and post-processing. 

The algorithmic steps 2-6 are fully automated, i.e. extensive user interaction is not required for 
these steps. In this section these algorithmic steps are discussed. 

2.1 Geometry definition 

The geometry of an aircraft configuration is usually defined using a CAD system (e.g. ICEM- 
CFD) and described in a standard data format. Typically, the geometry is established in terms of 
geometric entities, viz. support surfaces and curves. The relation between these geometric entities, 
referred to as the topology, has to be defined explicitly. The respective entities are then modelled 
by means of piecewise cubic surface and curve representations. 

At present an application starts either with a user-defined aircraft configuration or an aircraft 
configuration abstracted from a three-dimensional multi-block grid [4]. By adding the geometry 
description of far field boundaries and the symmetry plane (optional) a three-dimensional flow 
domain enclosed by bounding surfaces is defined. In such a way the boundaries of the flow domain 
are fully defined. 


2.2 Surface triangulation 

In order to perform a flow calculation the flow domain has to be diseretised, i.e. boundaries are 
triangulated and a tetrahedral grid is generated in the interior of the flow domain. Firstly, the 
boundaries of the flow domain are diseretised resulting in a surface triangulation. To this purpose 
a surface triangulation algorithm is employed which incorporates an equi-distribution algorithm for 
discretising curves and an advancing front-type generation algorithm for discretising surfaces. A 
distribution function controls the size of the edges and triangles in the surface triangulation. 

This distribution function is defined by parameters which are specified in the background grid and 
in source terms. Source terms are introduced in regions to refine regions of special interest. The 
distribution function expresses a desired spacing in each point of the flow domain. 

2.3 Three-dimensional grid generation 

Subsequently, a three-dimensional grid of Delaunay-type is obtained by employing a three-dinu nsional 
grid generation algorithm. In this algorithm a prospective node is located at the bary center of a 
tetrahedron. If this prospective node satisfies the Delaunay criterion the node is inserted and con- 
nected to the existing tetrahedral grid. In case the grid size, due to insertion of this node, becomes 
to small the prospective node is rejected. Smoothness of the three-dimensional grid is treated 
explicitly by the distribution function (as defined by the background grid and source terms). As 
a result a three-dimensional tetrahedral grid is generated which is bounded by a topologically 
two-dimensional surface triangulation. 


2.4 Pre-processing 
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The next algorithmic step is the pre-processing algorithm. This algorithm is designed to optimise 
the flow calculation for the generated three-dimensional grid. By employing a colouring algorithm 
for the edges in the three-dimensional grid a high degree of vectorisation in the flow calculation 
can be reached. 

2.5 Flow calculation 

In the flow calculation algorithm three-dimensional inviscid flow is simulated. 1 he incorporation of 
an upwind solver makes the CFD-system suited for the simulation of subsonic, transonic, supersonic 
as well as hypersonic flows. For the spatial discretisation of the three-dimensional Luler equations 
Roe's approximate Riemann solver [5] is utilised. Second-order accuracy is achieved by employing 
a MUSCL interpolation [6] for the state vectors, following the approach in [7, S], An entropy 
fix is incorporated to prevent physically incorrect features in a numerical solution. At present 
the boundary conditions are treated with first-order accuracy. Time integration is established by 
adopting a Runge-Kutta time step algorithm. Convergence acceleration is achieved bv local time 
stepping and residual averaging. 


2.6 Grid adaption 

In case the flow calculation has sufficiently converged (by taking a large number of Runge-Kutta 
time steps) the grid adaption algorithm can be employed. The grid adaption algorithm is based 
on remeshing. Based on the values of an adaption variable (for instance: density, Mach number or 
pressure coefficient) an adaption indicator is calculated which is based on the undivided difference of 
this adaption variable. According to the value of the adaption indicator sources terms are generated 
which are utilised to adapt the surface triangulation. Subsequently, the three-dimensional grid 
generation algorithm is adopted to generate an entire new three-dimensional grid. By incorporating 
the newly generated source terms in the existing distribution function a locally refined grid is 
obtained. For the adapted grid a numerical solution is then obtained by first employing the pre- 
processing algorithm and subsequently the flow calculation algorithm. This process can be repeated 
until a numerical solution of sufficient accuracy is obtained. 

2.7 Visualisation and post-processing 

Finally, a post -processing algorithm is adopted to calculate aerodynamic quantities. Flow visuali- 
sation is achieved by using the commercial package Data Visualiser [9], Interfaces are available to 
visualise the geometry definition, the surface triangulation and the numerical solution. 

3. APPLICATIONS 

This section discusses two applications of the current (T I) system at NLR. The first application 
concerns transonic flow around a wing/body configuration and the second application concerns 
hypersonic flow for a generic three-dimensional mixed-compression airbreathing inlet configuration. 

3.1 AS28G wing/body configuration 

In the first application described here inviscid flow around the AS28G wing-body configuration is 
simulated. The geometry is defined by 14 support surfaces and 42 curves. The physical coordinates 
of these support surfaces and curves are abstracted from a multi-block grid. The support surfaces 
contain approximately 17000 nodes. Far field boundaries and a symmetry plane are incorporated 
to define a closed flow domain. 

The surface triangulation shown in figure 2 is obtained. It can be observed that the nose region, 
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tail region, wing leading edge, trailing edge and the tip are refined. Source terms are located in 
these regions in order to obtain sufficient resolution of the geometry of the AS28G wing/body 
configuration and to reduce expected losses in total pressure. 

One successive refinement is employed. The surface triangulation of the adapted grid is shown in 
figure 3. The dimensions of the surface triangulation and the three-dimensional grid for the initial 
and the adapted grid can be found in table 1. The problem size for the adapted grid is approximately 
two times larger in comparison with the initial. Cpu-times for the respective algorithmic steps are 
shown in table 2 and the convergence history is depicted in figure 1. 

In the flow calculation a three-step Runge-Kutta time step scheme with a CFL-number of 0.9 is 
utilised. Two jacobi iterations are employed to smooth the residuals. The flow calculation algorithm 
uses 81 words per grid node, which is relatively large [10], and reaches a performance of 253 Mflops 
on the NEC SX3. 

The Mach-number and total pressure distribution are shown in figures 4 and 5 respectively. A 
relatively small amount of grid points is necessary to resolve the lambda shock on the upper side 
of the wing. This is a promising result as for instance in a multi-block structured code more 
grid points are necessary to capture the same shock structure. Nevertheless, large losses in total 
pressure (±25%) are experienced at the leading edge of the wing. This is mainly due to fact, that the 
boundary condition at the wing is only first-order accurate. A relative large loss in total pressure 
is also observed on the upper side of the wing (8%,). 

3.2 Three-dimensional mixed compression airbreathing inlet 

One of the most critical enabling technologies of advanced reusable launchers is the propulsion 
system. For example a two-stage-to-orbit aerospace vehicle needs a propulsion system to power the 
vehicle from take-off to sustained flight at Mach numbers ranging from 6 to 7 (separation Mach 
number of the rocket powered second stage from the airbreathing first stage). Sustained flight of 
airbreathing aerospace vehicles at these Mach numbers is not feasible with todays technologies. 
The inlet considered here is of mixed compression type [11, 12]. before entering the inlet, the 
oncoming air is decelerated to a hypersonic Mach number lower than the flight Mach number by 
oblique shock waves (see figure 6). Near t he narrowest duct cross section, the throat, the flow passes 
a normal shock wave. Behind this shock the flow is subsonic. Further deceleration is required to 
a velocity acceptable for the combustion chamber. The subsonic deceleration requires a diverging 
duct called the diffuser. 

The present section discusses typical results for flow through the mixed compression airbreathing 
inlet. The geometry of the inlet is described in full detail in reference [13]. The inlet geometry is 
based on a two-dimensional design that has been obtained by using engineering tools. The inlet 
considered here is designed for a freest ream Mach number of A/ <XJ = 4.5. The forebody of the space 
plane is not modelled. Hence, a Mach number smaller than the flight Mach number is considered. 
Bleed slots are not taken into account since viscosity is neglected. Both the cowl lip and the side 
w'all are modelled with finite dimensions. 

The geometry of the inlet is defined by 46 curves and 18 support surfaces. The geometry contains 
16 planes and 2 cubic polynomial representations, namely the cowl lip and the leading edge of the 
side wall. Figure 7 shows an inside view' of the inlet from the upstream direction. The cowl lip 
and the side walls are clearly visible. At the top of the figure two ramps can be observed. The 
large bounding planes at the right-hand and left-hand side are boundaries of the flow domain on 
which symmetry conditions are imposed. In the inlet duct the right-hand bounding plane is also 
defined as a symmetry plane. At the left-hand side the wall of the inlet duct can be observed. An 
impression of the surface triangulation on the side-wall (left-hand side of figure 7) is shown in figure 
8. Sufficient grid points are generated to obtain an accurate representation of the cowl lip. 
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For the freest ream Mach-number M <Xl — 4.5 a typical Tinstart” condition is obtained, see figure 9. 
Tliis means that there is a strong bow shock in front of the duct inflow plane providing subsonic 
flow throughout the duct. In windtunnel experiments a ’’started” condition is obtained by opening 
the throat significantly so that the shock can pass. Once the shock has been swallowed, the diffuser 
throat area is decreased, while simultaneously the pressure ratio of the pressure of the outlet of the 
diffuser with respect to the freest roam pressure is increased so that the shock is caused to move to 
the diffuser throat [14]. 

During numerical simulations, however, the variation of the geometry is not trivia] to simulate. To 
solve this problem the calculation is started with a higher Mach number A/ x , = 5.0 so that effective 
throat area is increased. This results in a globally correct pattern of oblique shock waves. The 
unstart condition is not experienced as the effective throat area is too wide for that particular 
Much number. Subsequently, the Mach number is gradually lowered to the required value A/ x = 
1.5. 

Moreover, to obtain the final shock in the diffuser throat and accordingly subsonic outflow, the 
pressure at 1 he diffuser outlet has to be imposed. The location of the final shock is a function of 
the value of the pressure, which can be found by employing normal shock relations and using the 
fact that, the flow is isentropic behind the final shock (see reference [15]). 

figure 10 shows aside view of the "started" condition with a final shock behind the throat of the 
dint. It (an be observed that the oblique shocks and the final shock are resolved quite accurately. 

1. LEVEL OF AUTOMATION 

A reduction in U1 I)-problem-t urnaround time can be achieved by an increased level of automation. 

I his section will focus on the automation level of the current (TD system for the applications 
dismissed in section 4. I hose algorithmic steps in the current (T I) system described in section 2 
which require user-interaction are identified. This identification is necessary in order to assess the 
potent ial for an increased level of automation. 

A With respect to geometry definition (step 1) an aircraft configuration is usually defined in 
terms of standard CAD data format. Extensive user interaction is required to generate a 
multi-block grid (such as for the wing/body configuration) which is suited as input for the 
current (TD-syst.em. ( Wnsequontly the problem-turnaround-time is large for multi-block 
based geometries. 

B I he definition of far field boundaries (step 1) (such as for the wing/body configuration ) 
requires user interaction. I he topology of a three-dimensional far field cube, the dimensions 
of the cube and the relation with between cube and the aircraft configuration have to be 
established . 

C 1 he definition of the distribution function (steps 2 and 5), which controls the size and shape of 
ti i angles in the* surface triangulation and tetrahedral elements in the three-dimensional grid, 
ipquires user interaction. Visualisation of the generated surface triangulation is necessary 
in order to inspect; the influence of the distribution function. A surface triangulation is 
accepted if sufficient resolution of details of the aircraft configuration (like wing leading and 
trading edge) is obtained. After modification of the distribution function it is necessary to 
regenerate the surface triangulation. It is observed that the definition of a suitable distribution 
function does not guarantee a successful three-dimensional grid generation algorithm. An 
intolerably large number of tetrahedral elements may be generated in case parameters defining 
the distribution function are chosen too small. For the current (FT)-system an expert user 
is required to specify the distribution function. 
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I) It is necessary to define boundary types (like symmetry plane, far field or solid wall) for 
each surface of the flow domain (step 4). This requires only a relative small amount of user 
interaction. 

E This requirement also holds for the specification of flow calculation parameters (step 5), such 
as for instance: Mach-number, angle of attack, (T L-number, number of Runge-Kutta, time 
steps and Restart-option. 

F Inspection of the flow solution by visualisation (step 7) can have a negative influence on 
the problem-turnaround-time if the number of nodes in the three-dimensional grid is large. 
This effect can be attributed to the fact that the handling of a large data-file in a multi-usei 
environment and the operations necessary for visualisation, such as zooming, pivoting and 
moving, are relatively slow. 

G In the grid adaption algorithm (step 6) user interaction is required to assess the level of grid 
refinement. Initially, default values for the grid adaption parameters are taken. Visual in- 
spection then learns whether the surface triangulation has been adapted sufficiently. If the 
generated surface triangulation is not acceptable (for instance the level of refinement is too 
high) the grid adaption algorithm has to be repeated with modified input-parameters. More- 
over, the number of nodes generated in the three-dimensional grid (after adaption, remeshing) 
is very sensitive to the values of the parameters specified in the grid adaption algorithm. If 
the number of nodes generated becomes too large it is necessary to repeat the grid adaption 
algorithm and the three-dimensional grid generation algorithm. 

In the near future the focus is on the respective parts requiring extensive user interaction, namely 
items A, C and G. It is foreseen that the remaining items B, D and F requiring more than minimal 
user interaction can be automated. 

By constructing an interface with a standard CAD data format the generation of a multi-block 
grid as input for the CFD-system can be avoided (item A). The data structure in the current CFD 
system is capable of handling a standard CAD data format. This would already significantly reduce 
1 he amount of user interaction. 

The influence of the distribution function (background grid and source terms) and the grid adaption 
parameters on the three-dimensional grid generation algorithm remains difficult to estimate (item 
C and G). Ideally, a functional relationship between these parameters and the number of nodes 
generated in the three-dimensional grid could contribute here. However, for complex geometries 
this relationship is difficult to assess. 


5. CONCLUSIONS 

The capabilities of the CFD system have been demonstrated for two applications. For the AS28C 
wing/body configuration and the three-dimensional mixed compression airbreathing inlet configu- 
ration the results so far indicate the potential for short turnaround times. An assessment of those 
parts in the CFD system which require user interaction shows that the automation level can be 
increased. The actual problem- turnaround time (geometry definition, surface triangulation, dl) 
grid generation, flow calculation, grid adaption, flow visualisation and post-processing) lies in the 
order of a few days for both applications. 

Moreover, the results indicate that the current CFD-system produces accurate solutions. The 
upwind-based flow calculation algorithm yields a spatially second-order accurate solution. How- 
ever, the boundary conditions are currently only first-order accurate. 
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Surface triang. 

3D Grid Generation 

Flow calc. 

Grid 

nodes 

triangles 

nodes 

elements 

edges 

time steps 

Initial 

22275 

44546 

105786 

609338 

737396 

500 

Adapted 

30744 

6 1484 

222359 

1348274 

1600926 

2000 


Table 1: Number of nodes and triangles in the surface triangulation, number of nodes and number 
of elements in three-dimensional grid for the initial grid and the adapted grid for the AS2XG 
wing/body configuration 


Algorithm 

cpu-time 
(initial grid) 

cpu-time 
(adapted grid) 

M flops 

computer 

Surface triangulation 

31m 20s 

- 

4 

SGI-Onyx 

3D grid generation 

5m 54 s 

14m 16s 

4 

NEC SX3 

Pre-processing 

lm 18s 

3 m 3 s 

3 

NEC SX3 

Flow calculation 

46m 11s 

8h 5m 6s 

253 

NEC SX3 

Grid adaption 

2m 58s 

- 

5 

NEC SX3 


Table 2: Cpu-times for the respective algorithmic steps for the AS28G wing/body configuration 



Figure 1: Convergence history of the flow calculation algorithm for the AS28G wing/body config- 
uration on the initial grid and the adapted grid 
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Figure A: Surface triangulation for Uie AS 28 CJ wing/body configuration after one successive refine- 
ment 
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Figure 4: Mach number distribution for the AS28G wing/body configuration on the upper side of 
the wing after one successive refinement 
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Figure 5: Total pressure distribution scaled with respect to the freestream values for the AS28G 
wing/body configuration on the upper side of the wing after one successive refinement 


724 




Figure (j: Main inlet components and terminology 




Figure 7: Impression of the three-dimensional mixed compression airbreathing inlet geometry mod- 
elled by cubic surface and curve representations 
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Figure 9: Side view of the Mach number distribution. The flow condition represents ’’unstarted” 
flow: strong bow shock and subsonic flow in duct 
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Figure 10: Side view of the Mach number distribution. The flow condition represents ’’started” 
flow with a final normal shock in the duct 
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Figure 11: Detailed view of the Mach-number distribution at the intersection of the side- wall ramp 
and the cowl lip. The flow condition represents "started flow” with a final shock in the duct 
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